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High-order methods are quickly becoming popular for turbulent flows as the amount 
of computer processing power increases. The flux reconstruction (FR) method presents a 
unifying framework for a wide class of high-order methods including discontinuous Galerkin 
(DG), Spectral Difference (SD), and Spectral Volume (SV). It offers a simple, efficient, and 
easy way to implement nodal-based methods that are derived via the differential form of the 
governing equations. Whereas high-order methods have enjoyed recent success, they have 
been known to introduce numerical instabilities due to polynomial aliasing when applied 
to under-resolved nonlinear problems. Aliasing errors have been extensively studied in 
reference to DG methods; however, their study regarding FR methods has mostly been 
limited to the selection of the nodal points used within each cell. Here, we extend some 
of the de-aliasing techniques used for DG methods, primarily over-integration, to the FR 
framework. Our results show that over-integration does remove aliasing errors but may 
not remove all instabilities caused by insufficient resolution (for FR as well as DG). 

I. Introduction 

Accurate and efficient numerical methods for the simulation of complex aerodynamic flows are desired 
to help reduce the design and development costs of aerodynamic vehicles and aeropropulsion systems. The 
primary challenge lies in the prediction of the turbulent motion within the flow and its effect on the overall 
motion of the fluid. Advanced methods such as large-eddy simulation (LES) and direct numerical simulation 
(DNS) offer great advancements in predicting turbulence, but these methods require significant improvements 
in accuracy and efficiency in order to become widely used. Focusing on either the accuracy or efficiency aspect 
of the underlying numerical algorithms, when developing a computational fluid dynamics (CFD) code, usually 
demands concessions from the other parameter. 

Promising results have been found applying LES in the realm of aeropropulsion, but the spatial and tem- 
poral resolution needed for LES requires high-order (higher than second) numerical schemes historically only 
found in structured grid methods with large stencils. These methods have numerous drawbacks: difficulty in 
grid generation, high sensitivity to grid quality, complex boundary conditions, and poor scalability for paral- 
lel computing. These drawbacks severely limit the ability of these codes to compute complex aeropropulsion 
configurations, e.g., noise suppressing nozzles . * 1 * 

In an attempt to alleviate the difficulties of structured grid generation for complex geometries, researchers 
have applied LES methods to unstructured grids . 2,3 However, the limited second-order accuracy of current 
unstructured methods requires an exorbitant number of grid points and offers little incentive over the tra- 
ditional structured methods . 4,0 The flux reconstruction (FR) method 6 * 8 provides a simple, efficient, and 
easy to implement method for solving the Navier-Stokes equations on unstructured grids and is accurate to 
an arbitrary order. By employing the high-order FR method with LES on unstructured grids, limitations 
related to structured grid generation, grid resolution, and computational efficiency can be overcome . 9 * 12 

High-order spectral- /finite-element type methods, including discontinuous Galerkin (DG) and FR meth- 
ods, use a finite set of basis functions to construct a numerical solution within each grid cell. When applied 
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to nonlinear equations, all of these methods can generate errors if any of the modes arising from the nonlin- 
ear terms are outside the span of the set of basis functions. In these cases, the energy from the unresolved 
(under-integrated) modes are aliased onto the lower modes producing so-called aliasing errors. The effect 
aliasing errors can have on a given simulation is highly unpredictable, and many researchers have shown 
they can produce severe instabilities within a problem. 13-21 

Aliasing errors have been extensively studied in reference to DG methods resulting in numerous stabiliza- 
tion techniques aimed at de-aliasing the numerical solution. One of the most commonly used techniques is 
a modal filter 16,18 ’ 22 which is used to add dissipation to the highest modes, thus reducing the magnitude of 
the oscillations within the polynomial. However, the common exponential filter contains three adjustable pa- 
rameters and the amount of filtering required to remove any instabilities is highly dependent on the problem 
being solved. Therefore, filtering in general is not very user friendly. Other stabilization methods include but 
are not limited to: the spectral vanishing viscosity method 13,19 ’ 23-2 ' (artificial viscosity), superconsistent 
collocation, 28,29 and polynomial de-aliasing 13,14 ’ 18 ’ 19 which is also referred to as over-integration. 

Regarding FR methods, the study of aliasing errors has mostly been limited to the selection of the nodal 
points used within each cell. 17 ’ 20 ’ 21,30-33 In this work, we look to extend some of the de-aliasing techniques 
used for DG methods, primarily over-integration, to the FR framework. In section II, we briefly review the 
FR method for the ID conservation law. In section III, we derive the interpolation and projection operators 
and subsequently describe their use for over-integration within the FR method. Section IV describes the 
governing Euler equations. In section V, we apply over-integration to many different cases of the isentropic 
Euler vortex problem to see the effects it has on simulations of varying resolutions. Finally, conclusions and 
discussion are provided in section VI. 


II. Numerical Methods 


II. A. Spatial Discretization 

We briefly review the FR method below. Consider the ID conservation law 

Ut fx = 0 (1) 

where the subscripts t and x denote partial differentiation with respect to time and space, respectively. 
Let the computational domain fl be divided into an arbitrary number of nonoverlapping cells. Each cell is 
transformed from physical space within the global domain to a common reference space/element, I, which 
in ID is the line interval, I = [—1, 1]. With r denoting the local coordinate variable in the reference space, 
the mapping function r = x(x) is used to complete this transformation. Similarly, the inverse function 
x = X _1 (r - ) maps the reference element back to physical space. The details of such mappings, especially for 
higher dimensions, are beyond the scope of this paper; the interested reader can find additional information 
in many texts on spectral-/finite-element type methods, including but not limited to Ref. 15, 16, 19,34. The 
following review focuses on using the FR method to solve the transformed version of equation (1) within 
each reference element. 

At a given time t , let the exact solution to the transformed conservation law within the reference element 
be given by u E (r,t ) which may be a nonsmooth and/or discontinuous function. Let u(r,t ) be a polynomial 
of degree P and approximate u E (r,t) within the reference element. This approximate solution polynomial , 
rt(r, t), is local to each element, and is generally discontinuous across cell interfaces. With the dependence 
on time being understood, we will drop the variable t from u E {r,t) and u(r,t) to simplify notation. 

Let Np refer to the minimum number of solution points needed to define the solution polynomial within 
the reference element, and let £ denote the set of solution points, n = 1,2,..., Np. Using tensor products 
to extend this to higher dimensions, the jV-dimensional conservation law requires Np = (P + 1)^ solution 
points within each computational cell. 

The Lagrange polynomial is defined as the interpolating polynomial through a set of nodal points, for 
example that takes the value 1 at £ n and zero at all other points 

n p _ p 

Ur ) = n f -= T ( 2 ) 

i= 1 ^ 
i^n 

Note that the Lagrange polynomial in equation (2) is a polynomial of degree P = Np — 1 and the set of 
Lagrange polynomials for all solution points in £ is a set of Np polynomials each of degree P. 
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A nodal coefficient of the solution, denoted by u n , is simply the value of the approximate solution 
polynomial evaluated at the solution point, i.e., u n — u(f n ). Let u P denote a vector containing the set 
of Np nodal coefficients for the solution at each of the solution points. Assuming that all of u P is known, 
we can use the Kronecker delta property of equation (2) to construct the approximate solution polynomial 
as 

N P 

u(r ) = y' j u n i n (r) (3) 

n—1 

The above relation is why the set of Lagrange polynomials can also be referred to as a nodal basis set. In 
this section and section III, the tilde accent will be used exclusively to identify the nodal coefficients for a 
polynomial function. 

The derivative of the solution polynomial is a linear combination of the derivatives of the polynomial 
basis functions 

d Np 

— u(r) = u r (r) = y^u n 

n—1 

Using the previous equation, we can define the derivative operator, D r , as the Np x Np matrix with elements 


dr 


( 4 ) 


D r [ i,j] 



( 5 ) 


that transforms the nodal coefficients, u P , into the derivative nodal coefficients, u r . Note that this derivative 
is exact only when the function being differentiated is a polynomial of degree P or less. 

If we use the nodal solution to evaluate the flux function, /, at each solution point, then interpolate 
these Np values by a polynomial of degree P, we get what is called the discontinuous flux function. It only 
uses information inside a given cell and does not account for interaction between it and the data in the 
neighboring cells. Using the superscript ‘D’ for discontinuous, the discontinuous flux polynomial, / D (r), is 
given by 

N P 

/ D ( r ) = 5Z f^nfr) (6) 

n—1 

where the nodal coefficients of the discontinuous flux, f n , are evaluated at the solution points using the 
nodal coefficients of the solution, i.e., /„ = f(u n ). This collocation projection results in / D (r) being a 
polynomial of degree P, which only matches the true flux function if the conservation law in equation (1) is 
the linear advection equation. For a nonlinear conservation law such as those found in the Euler equations 
or Navier-Stokes equations, the true flux function evaluated using the solution polynomial is of much higher 
degree than P and possibly even irrational, and / D (r) above may not provide a good approximation. We 
will explore other ways of approximating the true flux function later. 

With / D (r) given by equation (6), we can compute the derivative of the approximated discontinuous flux 
using the derivative operator from equation (5). Denoting / as the column vector containing the coefficients 
/„, n = 1,2,..., Np, the derivative ^ [/ D (r)] can be approximated using the matrix- vector product 


fr = Drf (7) 

where, f r is the column vector containing the values of the approximate flux derivative at each solution point. 
Note that the collocation projection of the discontinuous flux in equation (6) is of degree P. Consequently, 
its derivative is of degree P — 1. 

The method for computing the divergence of the discontinuous flux given by equation (7) will be referred 
to as the Lagrange polynomial (LP) method since it utilizes the derivative of the Lagrange polynomials to 
evaluate ff. An alternative method can be formulated by using the chain rule to rewrite equation (7) as 

fr = | ^D r ii P ( 8 ) 

ou 

Naturally, this will be referred to as the chain rule (CR) method and requires the evaluation of the flux 
Jacobian at each solution point, making it more expensive than the LP method. Whereas numerical exper- 
iments have shown that the CR method yields more accurate solutions over the LP method for nonlinear 
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problems, 3 ' 1 the method itself is not fully conservative. A modification to the CR method was presented in 
Ref. 36 that ensures conservation. 

If we employ either of the above methods for ff to evaluate f x for the conservation law, we obtain 
erroneous solutions since such a derivative does not include the interaction of the data between adjacent 
cells. Therefore, we seek to construct a so-called continuous flux function, denoted by F(r), which accounts 
for this interaction between adjacent cells and approximates the discontinuous flux function in some sense, 
to which we can then calculate its derivative. The continuous flux function will be obtained by adding a 
correction to / D , the discontinuous one. 

To this end, at each (cell) interface, let u L and u R be the values taken on by the solution polynomials 
to its left and right, respectively. An upwind flux value common for the two adjacent cells, denoted by 
/* = /*(ul,Ur), can then be defined via the wind direction. Here, for the case of the Euler equations, we 
use Roe’s flux 3 ' with an entropy fix. 38 

Let f BB denote the upwind flux at the ‘left boundary’ of the cell and /* B denote the upwind flux at the 
‘right boundary’ of the cell. We can then construct the continuous flux function, F(r), as a polynomial of 
degree P+1 by requiring that it takes on the upwind flux values at the two cell interfaces, and it approximates 
/ D via P — 1 conditions. Instead of defining P, we define F — f D , which approximates the zero function. At 
the two interfaces r = ±1, the function F(r) — f D {r) takes on the following values, which are called the flux 
jumps: 

p(-i) - n- 1) = f* b - a- 1) and p(d - ni) = /* B - m (9) 

Therefore, F(r) — / D (r) has prescribed left and right interface values, is of degree P + 1, and approximates 
zero. 

We now separate the prescription of the jump at the left interface from that of the right. This separation 
plays a critical role in the new approach. Again, using the local coordinate, let g LB be the correction function 
for the left interface of the reference element defined by 

g LB (— 1) = 1, <7 lb(1) = 0 (10) 

and g LB is a polynomial of degree P + 1 approximating the zero function in some sense. We always choose 
g RB by reflection: g RB (r) = g LB (— r). As a result, 

<7rb( — 1) = 0, Srb(1) = 1 (11) 

Consider the left interface r = — 1. The polynomial 

/ D (r) + [/*B-/ D (-l)]ffBB(r) (12) 

provides a correction at the left interface for / D (r) by changing the flux value at this interface from / D (— 1) 
to f* g, while leaving the value at the right interface unchanged, namely, / D (l). Next, the polynomial 

F(r) = nr) + [/lb - / D (-l)] dUr) + [/ RB - / D (l)] SrbM (13) 

provides corrections to both interfaces; using equations (10) and (11), one can easily verify that F(— 1) = f BB 
and F(l) = / RB . Thus, the above F(r ) is of degree P + 1, takes on the upwind flux values at the two 
interfaces, and approximates / D (r) in the same sense that g L B and g RB approximate the zero function. 

The derivative of F(r) at the solution point £„ is evaluated as 

= F r (Z n ) = f r , n + [f* B - / D (-l)] <?' B (U + [/r B - / D (l)] ^ B (U (14) 

U 

Finally, F r is transformed from reference space to physical space using the inverse mapping X -1 , where it 
can then be employed to approximate f x at each solution point within the physical grid cell. The solution 
to equation (1) is then advanced in time via an explicit Runge-Kutta method. 

In order to complete this brief review of the FR method, we must define the correction function g BB 5 the 
definition of g RB follows by reflection as previously discussed. With L k denoting the (standard) Legendre 
polynomial of degree k, the right Radau polynomial of degree k (with k > 0) is defined by 

RR,k = ( ^ ( L k — pfc-i) (15) 
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Figure 1: An example of polynomial aliasing when computing a nonlinear flux function. The blue curve is 
the solution polynomial and the blue dots give the value of this polynomial at the solution points. The red 
curve is the analytic flux polynomial and the green marks give the flux function evaluated at the solution 
points. A collocation projection results in the green curve which is the “aliased” flux polynomial through 
the flux values at the solution points. The brown curve is the L 2 projection of the analytic flux onto the 
solution polynomial space. 


A family of correction functions for the left boundary that results in stable schemes is given by 

g LB = a R.R'P +1 + (1 - a) R R: p with 0 < a < 1 (16) 

Note that a = 1 results in the DG method, and a = 2 P +1 results in the modified Spectral Difference (SD) 
method that is stable. 

II. B. Temporal Discretization 

The results presented in this paper were integrated in time using the 3-stage strong stability preserving (SSP) 
Runge-Kutta method, 39 which is 3rd-order accurate in time. All of the results used a constant (global) time 
step to advance the solution in time. The size of the time step for each simulation was chosen to be sufficiently 
small meaning further reduction in its size would not produce significant changes in the solution error. 

III. De-Aliasing using Over-Integration 

High-order spectrah/finite-element type methods, including DG and FR, use a finite set of basis functions 
to construct a numerical solution within each grid cell. When applied to nonlinear equations, all of these 
methods can generate errors if any of the modes arising from the nonlinear terms are outside the span of 
the set of basis functions. In these cases, the energy from the unresolved modes are aliased onto the lower 
modes producing so-called aliasing errors. The effect these errors can have on a given simulation is highly 
unpredictable, and many researchers have shown they can produce severe instabilities within a problem. 13-21 

An example of how these aliasing errors can develop is shown in figure 1 where we have defined the flux 
function to be the square of the solution, i.e. , f(u) = u 2 . The blue curve represents the solution within 
the grid cell as a polynomial of degree 4 and the five blue dots give the value of the solution at each of the 
solution points. Analytically applying the flux function to the solution polynomial results in the analytic 
flux polynomial, which is shown as the red curve and is a polynomial of degree 8. Evaluating the analytic 
flux polynomial at each of the solution points results in the green ‘x’ symbols. A collocation projection of 
the analytic flux is the polynomial of degree 4 that interpolates these values of the flux function evaluated 
at the solution points. This is represented as the green curve in figure 1 and illustrates how the analytic 
flux polynomial is “aliased” onto a polynomial space of lower degree. Finally, a Galerkin or L 2 projection of 
the analytic flux onto the solution polynomial space is the best possible approximation in the L 2 norm that 
we can achieve of the analytic flux as a polynomial of degree 4. This projected flux polynomial is given by 
the brown curve and the brown squares mark the value of the polynomial at the solution points in order to 
show the difference from the collocation flux values. 

One way to prevent aliasing errors is through over-integration, which is also referred to as polynomial 
de-aliasing. Let the analytic flux , / A , be the exact analytic expression resulting from the evaluation of the 
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flux function using the solution polynomial u , i.e. , f A = f{u). The goal of over-integration is to get an 
optimal approximation of the analytic flux as a polynomial within the solution polynomial space. This is 
accomplished by using a set of quadrature points to increase the sampling of the flux function above what 
is capable by the solution points. The interpolating polynomial created from this increased sampling is of a 
higher degree than the solution polynomial space. Finally, applying an L 2 projection to this higher-degree 
flux polynomial results in the best possible approximation of analytic flux in the L 2 norm as a polynomial 
within the solution polynomial space. 

If the analytic flux happens to be a polynomial, we can recreate it with minimal error if enough sampling 
points are used. Figure 1 is an example of the analytic flux being a polynomial. For this particular example, 
if we used the value of the solution polynomial at 9 optimally-located points to evaluate the flux function, 
we could recreate the red curve exactly. However, care must be taken to limit the amount of over-integration 
that is used to only what is necessary to exactly recreate the analytic flux. At some limit, further sampling 
becomes superfluous, wasting computational resources to yield no additional information. 

On the other hand, if the analytic flux is not a smooth and/or continuous function, no amount of over- 
sampling will be able to exactly recreate the analytic flux as a polynomial, and some aliasing error will be 
introduced because of this approximation. In these situations, we choose an amount of over-integration and 
hopefully the aliasing error will be negligible compared to the overall simulation error. 

We implemented over-integration into our FR solver in two different ways. The first method uses the 
quadrature points to over-sample the flux function resulting in an approximation of the analytic flux as a 
polynomial of degree greater than P. An L 2 projection is then used to project this over-sampled flux onto 
the solution polynomial space. This method will be referred to as flux over-integration, and it involves a 
complex series of steps that adds two interpolation and two projection operations to the base FR method. 
The second method uses the quadrature points to over-sample the residual function to get an optimal 
approximation of the residual within the solution polynomial space. Referred to as residual over-integration, 
this method is much easier to implement since it adds only one projection operation to the base FR method. 
However, residual over-integration is generally more expensive since it requires computing the complete 
residual function using the quadrature points. Further details on both of these methods are provided in 
sections III.B and III.C. 

III. A. Interpolation and Projection Operators 

In order to use over-integration, we need to define two matrix operators: the interpolation and projection 
operators. The projection operator, P n , is used to project a nodal solution onto the space of polynomials of 
degree n, i.e., T n . The interpolation operator, X v , is used to evaluate a function at a set of arbitrary nodal 
points, r). For the sake of brevity, some details will be omitted in the process of defining these two operators 
that is to follow. A more thorough explanation of the general procedure being presented here can be found 
in Ref. 19. 

We start by identifying some fundamental parameters and their respective notation. Recalling that P is 
the degree of the solution polynomial, let Q refer to the degree of the polynomial used for over-integration, 
where Q > P. Let Nq refer to the minimum number of quadrature points needed to define the over- 
integration polynomial within the reference element I. Using tensor products to extend the ID reference 
element to higher dimensions, the total number of quadrature points for the Af-dimensional reference element 
is Nq = (Q + 1)^. 

Note that we will use the term nodal , as in nodal set or nodal points, to refer to any arbitrary set of points 
within the reference element, e.g., both the solution and quadrature points are sets of nodal points. In order 
to discern between these two nodal sets we will adhere to the following notation. We retain the notation 
from section II. A when referencing the solution points: £ will denote the set of Np solution points, and 
the subscript P will denote a set relating to the solution points with the index n referring to a component 
of such a set, e.g., £ n . We define similar notation for the quadrature points: r) will denote the set of Nq 
quadrature points, and the subscript Q will denote a set relating to the quadrature points with the index 
q referring to a component of such a set, e.g., r] q . The term £ will be used to denote an arbitrary nodal 
set containing N nodal points, Q, whenever we need to present a concept that is not unique to either the 
solution or quadrature points. Finally, the subscripts i and j will represent dummy indices which will be 
used whenever one of the dedicated indices is already present within an equation. 

To help illustrate the various nodal sets used within the reference element, an example of a 2D reference 
element is provided in figure 2. This example corresponds to an element with P = 1 and Q = 4 resulting 
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Figure 2: A generic grid cell illustrating the solution 
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and 


in Np = (P + l ) 2 = 4 solution points (shown as blue squares) and Nq = {Q + l ) 2 = 25 quadrature points 
(shown as red ‘x’ symbols). Shown along the edges of the element in figure 2 are the interface/flux points 
required to account for interaction between adjacent cells. The nodal points along each edge correspond to 
the respective ID nodal sets, with the solution flux points shown as green squares and the quadrature flux 
points shown as black circles. 

The choice of nodal sets that are used can have an affect on the stability, accuracy, and efficiency of nodal- 
based polynomial methods like FR 31 and DG. 40,41 There are two families of nodal points that are used within 
this work: the Gauss-Legendre and the Lobatto-Gauss-Legendre nodes. In ID, the Gauss-Legendre nodes, 
or more simply the Gauss nodes, correspond to the N zeros of the Legendre polynomial of degree N within 
the interval [—1,1]. The Lobatto-Gauss-Legendre nodes, or more simply the Lobatto nodes, correspond to 
the N — 2 zeros of the derivative of the Legendre polynomial of degree N — 1 within the interval [—1,1], plus 
the two end points ±1. 

Each of these nodal sets is associated with a quadrature rule that defines a matching set of quadrature 
weights for each nodal point. Note that the terms quadrature rule and quadrature weights imply no direct 
reference to the quadrature points and simply refers to a set of nodes and matching weights that can be used 
to numerically approximate an exact integral by finite summation. Let w denote the set of N quadrature 
weights, Wi, i = 1,2, ...,1V, associated with each of the nodal points in £. With a(r) representing some 
arbitrary function defined on I, we can use the quadrature rule associated with £ to approximate the following 
integral as 

r i n 

/ a(r) dr Ki^2a(Ci)wi (17) 

J - 1 i=i 


The benefit of using the set of N Gauss nodes and associated quadrature weights is the approximation in 
equation (17) is exact when air) is a polynomial of degree 2N — 1 or less. When using the set of N Lobatto 
nodes, the approximation in equation (17) is exact when a(r) is a polynomial of degree 2 N — 3 or less. 

Next, the inner product between two arbitrary functions, a{r) and b{r), within the reference element I, 
is defined by the equation 


(a, 6) 


J a{r) b(r) dr 


(18) 


The two functions a(r ) and b{r) are considered orthogonal to each other if the integral in equation (18) 
evaluates to 0. Using the quadrature rule associated with £, equation (18) can be approximated as 


N 

(a, 5) ~ '52a(Q)biQ)w i 




(19) 


Let a and b denote vectors that contain the values of their respective functions, a (r) and &(?'), evaluated at 
each of the nodal points in £. Additionally, let w ^ denote the set of quadrature weights associated with £. 
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The weight matrix, W , is then defined as the square, diagonal matrix with w £ located along its diagonal, 
i.e. , W = diag(io^). Using a, 6, and W , the summation in equation (19) can be rewritten in matrix form as 

(a,b)=a T Wb (20) 


In equation (20), the assumptions have been made that a and b are both polynomials and the quadrature 
rule associated with the nodal set £ is sufficiently accurate to exactly evaluate the integral in equation (19). 

We wish to define an interpolation operator, X, which outputs the nodal coefficients at the quadrature 
points when given the nodal coefficients at the solution points. Let u q and u Q relate to the quadrature points 
and be the respective equivalents to u n and u P . The nodal coefficients for all the quadrature points, can be 
found using equation (3) 

N P 

u q = u(rj q ) = y ^u n t„{ri q ) for q = 1, 2, . . . , Nq (21) 

n—1 


It is clear that the Lagrange polynomial term on the right side of equation (21) is simply a Nq x Np dimen- 
sioned matrix that is multiplied by the vector u P resulting in u Q . Therefore, we can rewrite equation (21) 
in matrix form as 

Uq =■ X V U P (2*2) 


where the interpolation operator, X rj . with elements X qn = £ n (jl q ), is a matrix that interpolates the solution 
nodal coefficients to the quadrature points, rj. For example, if P = 1 and Q = 2, the interpolation operator 
is the matrix 


T — 

-*■'77 — 


kiivi) 

kAm) 

kAm) 


k 2 (vi) 
k 2 ) 
kAm). 


(23) 


where is used to specify the Lagrange polynomial that takes the value 1 at the solution point and zero 
at all other solution points k^n- 

Instead of using the nodal form of the solution polynomial given by equation (3) , we can decompose the 
polynomial into hierarchical modes by using a set of modal basis functions and their corresponding modal 
coefficients. For an aaD solution, let Nm refer to the maximum (or final) mode which corresponds to a 
polynomial basis function of degree a rxP. Whereas Nm denotes the maximum mode, the total number of 
modes is actually the same as the total number of solution points because we have defined the modal index, 
m, to begin at 0; this was done to signify that the first mode, i.e., m = 0, corresponds to the constant 
mode where all the polynomial exponents are 0. Thus, the value of Nm is Nm = (P + 1) A — 1 = Np — 1. 
Combining all of these notations, let u P denote the set of modal coefficients, u m , m = 0, 1, ... , Nm ; we have 
used the ‘hat’ accent in order to distinguish it from the tilde accented nodal coefficients. 

Let 4> P be the set of Np modal basis functions, </> m (r), to = 0, 1, ... , Nm, that directly correspond to the 
modal coefficients in u P . Note that for the ID case, the basis function </> m (r) is a polynomial of degree to. 
Let B denote the basis matrix where each entry in the matrix is given by B,j = Thus, each column 

vector of the matrix represents one of the basis functions evaluated at each point in a nodal set, whereas 
each row vector of the matrix represents each basis function evaluated at a specific nodal point. 

In the formulation of the projection operator, we will deal with two basis matrices. The first basis matrix, 
which retains the notation B , or has the dimensions Nq x Np and corresponds to each of the Np 

basis functions evaluated at each of the Nq quadrature points. Using the same example of P = 1 and Q = 2 
from equation (23), the basis matrix expands to 


B 2.1 = 


<M*7i) 

<M??2) 


0i(*7i) 

<M*72) 

0i (m). 


(24) 


In the present work, the normalized Legendre polynomials are used as the modal basis functions. These are 
ideal because they fit the hierarchical requirement and the Legendre polynomial of degree to is orthogonal 
fO *km— 1* 

The second basis matrix that we will use is denoted by V, and represents the Np basis functions evaluated 
at each of the Np solution points, i.e., V = Bp p. The matrix V is denoted differently than B to immediately 
recognize that it is square and invertible, whereas B is generally not square and thus not invertible. 
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Using the modal coefficients and corresponding basis functions, the approximate solution polynomial 
from equation (3) can be reconstructed using 

N m 

u{r ) = ^2 Um^mir) (25) 

m — 0 

If we use equation (25) to evaluate the nodal coefficients for each of the Nq quadrature points, we end up 
with 

Nm 

U q = u(r] q ) = ^2 UmfimiVq) for q = 1,2,..., Nq (26) 

m = 0 

Similar to how equation (21) was rewritten as equation (22), we can rewrite equation (26) in matrix form as 

u Q = Bu P (27) 

If we were instead to use equation (25) to evaluate the nodal coefficients for each of the Np solution points, 
we end up with the matrix equation 

Up = Vu P (28) 

Because the matrix V is invertible, we can very easily rearrange the above equation to solve for the modal 
coefficients given the nodal coefficients at the solution points, i.e. , 

up = V _1 u P (29) 

From the previous equations, we can see that the basis matrix is essentially an operator that transforms the 
modal solution into the nodal solution. This completes the definition of all the necessary components that 
are required for us to define the projection operator. 

Before we continue on to the projection operator, however, we take the opportunity to quickly present an 
alternative method for evaluating the interpolation operator from the basis matrices. If we plug equation (29) 
into equation (27), and then set the result equal to equation (22), we end up with 

B [V _1 m p ] = u Q = I v u P (30) 

The alternative formulation of the interpolation operator is clearly 

Z„ = BV~ 1 (31) 

We now continue on with defining the projection operator. Because the exact solution, u E (r ), is not 
necessarily in the space of polynomials spanned by the basis functions of the approximate solution, u(r), the 
residual error, R(u), from the approximation can be expressed as 

R{u) = u E (r) - u(r) (32) 

Let (p represent an arbitrary test function. We can apply the inner product of equation (32) with respect to 
the test function, <p, resulting in 

(<p,R(u)) = (<p,u E ) - (ip,u) (33) 

From the method of weighted residuals, we require that the left hand side of the previous equation is 0. 

The Galerkin projection, also referred to as the L 2 projection, uses the modal basis functions, 4 > , as the 
test functions in equation (33). Let us assume that we have used a set of quadrature points, rj, that is large 
enough to sufficiently sample the exact solution such that 

Nq 

u E (r)^^2u q e q {r) (34) 

9=1 

After setting the left hand side of equation (33) to zero, if we replace u E with the approximation in equa- 
tion (34), u with the modal formulation from equation (25), and the test function with the set of modal 
basis functions, we end up with 

Nm \ f Nq \ 

<t>i , T, ) = I Uq l q 1 for i = 0, 1, . . . , N M (35) 

m= 0 / \ 9=1 J 
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The summations can be pulled out of the inner products since u q and u m are just coefficients within their 
respective summations. Thus equation (35) simplifies to 

Nm Nq 

^ ^ Um [4*1 > 4*rn) — ^ ^ Uq [4'i •> ^g) for ^ — 0, 1 , . . . , Nm (36) 

m = 0 q = 1 

Using equation (19) and the quadrature rule associated with the quadrature points, the inner products in 
equation (36) can also be expressed as summations 

N m n Q n Q n Q 

^2 & (rij) fimivj) Wj = ^w g ^0i(?bKg(?7j)wb for i = 0, (37) 

m=0 j = 1 q= 1 j = 1 

Using equation (20), this can be rewritten in matrix form as 

B J WB u P = B J WI u Q (38) 

where I is the identity matrix created from the Lagrange polynomials for the quadrature points, i.e. , 


I = lq (rij) = Sqj for q, j = 1, 2, • . . , Nq (39) 

The leading term on the left hand side of equation (38), B J WB, results in a square matrix that can be 
inverted and moved to the right side giving 


Up = 


(b t WB^) 1 B J W 


Ur 


(40) 


The previous equation takes the set of Nq nodal values defined at the quadrature points, ii Q , and projects 
these values onto the solution polynomial space of degree P, resulting in the set of modal coefficients for the 
basis functions of the solution polynomial, u P . 

Using equation (28), we can transform the projected modal coefficients from equation (40) into the 
physical values at the solution points. Multiplying V by the term in square brackets from equation (40) 
gives the final projection operator, 

P p = v ( B J WB ) 1 B J W (41) 

where P P denotes that the projection is onto the polynomial space of degree P. Thus, the projection 
operator is used to project the Q-degree polynomial defined by the nodal coefficients u Q to a polynomial of 
degree P, which is then evaluated at the solution points resulting in the nodal coefficients u P , i.e., 


U P — P p Uq 


(42) 


Note that the projection operator above is simply a matrix and its use only requires multiplying it by the 
vector of nodal coefficients for the quadrature points. Since the matrices defining the operator are based on 
quadrature rules and basis functions that do not change during a simulation, the operator matrix is only 
required to be defined once during the preprocessing phase of the simulation. 

As a quick exercise, we take a look at what happens to the projection operator if over-integration is not 
used, i.e., Q = P and rj = £. In this specific case, the basis matrices B and V are equivalent. Substituting 
V for B throughout equation (41) results in 

P P = V (v T WV^j 1 V T W (43) 

Since all of these matrices are now square and invertible, we can use some matrix identities to expand 
equation (43) as 

P p = VV~ 1 W~ 1 (V T y 1 V T W (44) 

After canceling terms, the above projection operator is clearly the identity matrix. This is exactly what we 
expect because here we are trying to project a polynomial of degree P onto the polynomial space of degree 
P. 
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III.B. Flux Over-Integration 

The first over-integration method focuses on preventing aliasing errors that might arise from the evaluation 
of two flux related terms: the discontinuous flux divergence within the cell interior and the common flux 
at the cell interfaces. Because these two terms are independent of each other, we can split this method 
into two phases with each phase corresponding to one term. Both phases are similar in that they involve 
four basic steps: interpolating the nodal solution to the quadrature points, evaluate a flux function at the 
quadrature points, projecting this higher-order flux polynomial down to the solution polynomial space, and 
finally interpolating this projected flux polynomial from the quadrature points back to the solution points. 
As we showed in equations (41) and (42), the final step of interpolating the projected flux from the quadrature 
points to the solution points can be absorbed into the projection operator for improved efficiency since both 
operations are matrix multiplications. 

Let £ and r) respectively refer to sets of solution and quadrature points as in section III. A, where we were 
exclusively dealing with interior nodal points. However, we are now required to differentiate between interior 
and interface nodal sets, because the first phase of flux over-integration involves the flux function within 
the interior of the element, whereas the second phase involves the flux along the interfaces. We therefore 
introduce the subscripts ‘I’ and ‘F’ to denote the respective interior and interface nodal sets. The legend in 
figure 2 provides an example of this notation as applied to the reference quadrilateral element. 

The interpolation and projection operators for the interior points follow directly from section III. A. It is 
not immediately clear though, how these operators are applied to the interface nodal points where the spatial 
dimension is one lower than that of the interior nodal points. The interpolation operator is actually straight- 
forward; assuming we know the location of the interface points, the matrix that interpolates a polynomial 
from to either £ F or r] F is evaluated using a method similar to equation (21). Whereas the projection 
operator for the interface polynomial is still defined using equation (41), it requires all of the matrices on the 
right side of the equation to be recreated for the face reference element. In equation (41), the matrices V, 
B , and W are all defined for a a/-D element and not applicable to the (a/ — 1)-D polynomial we are looking 
to project at the interface. By using £ F and rj F to create interface versions of V, B , and W, we can create 
the interface projection operator with equation (41). 

III.B. 1. Over-Integrating the Discontinuous Flux Divergence 

The first phase of flux de-aliasing is applying over-integration to the evaluation of the discontinuous flux 
divergence. We start with assuming that the nodal coefficients are known at the interior solution points, £ : . 
The steps for evaluating the de-aliased discontinuous flux divergence are as follows: 

1. Interpolate the solution from to r) r 

2. Evaluate the flux function at r/ l using the interpolated solutions. 

3. Project the flux polynomial to the solution polynomial space. 

4. Interpolate the projected flux polynomial from r) l to £ : . 

5. Compute the derivative of the projected flux polynomial at each £ : . 

Following the previous steps results in the discontinuous flux as a polynomial of degree P , thus its divergence 
is of degree P — 1. This is same as the base FR method except here / D is a better approximation of the 
analytic flux in the L 2 norm. 

Additionally, we can define an alternative method to over-integrate the discontinuous flux by rearranging 
the final three steps as: 

3. Compute the derivative of the flux polynomial at each rj v 

4. Project the derivative of the flux polynomial to the solution polynomial space. 

5. Interpolate the projected flux derivative from r/ 1 to £ : . 

This slight modification results in /°(r) being a polynomial of degree P instead of P — 1 as it was before. 
However, this can create some consistency errors when evaluating the flux jumps at the cell interfaces. The 
values / D (±l), which are used to compute the flux jumps at the interfaces in equation (9), need to be 
consistent with the projected derivative. This means we need to create a discontinuous flux polynomial of 
degree P + 1 from f°(r) in order to compute the flux jumps in a consistent manner. 

Another consequence of projecting the derivative is that the derivative operator needs to be created from 
the Lagrange polynomials for the quadrature points instead of the solution points. Therefore, computing the 
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flux derivative at the quadrature points is significantly more expensive because of the increase in operations 
required by the matrix- vector product. Brainstorming these and other complications that we encountered 
led to an alternative method, residual over-integration, which is described in section III.C. 

III.B.2. Over-Integrating the Upwind Flux 

At each interface/flux point, the upwind flux, /*, is found by solving the Riemann problem that exists 
from the (generally) discontinuous solution between the two cells on the interface, i.e. , ul ^ ur- Numerous 
approximate Riemann solvers have been derived for solving this problem; they all operate on a pair of 
individual states and the underlying algorithms generally contain nonsmooth functions and/or conditional 
statements. However, for polynomial methods, the Riemann problem at an interface is actually between 
two discontinuous solution polynomials , tti(r) ^ UR(r). The solution to this Riemann problem is extremely 
complicated and all but guaranteed to not be a polynomial. Therefore, the upwind flux could very likely be 
the source of aliasing errors if we take an insufficient sampling of the total interface flux. 

The quadrature flux points, r/ F , are used to get a better sampling of the total interface flux, which can 
then be projected onto the solution polynomial space at the solution flux points £ p . For a given interface 
between two cells, the steps to achieve this are as follows: 

1. From the left and right cells on the interface, interpolate the solution from their respective to the 
shared quadrature flux points, r/ F . 

2. Evaluate /* at r) F using the left and right interpolated solutions. 

3. Project /* to the solution polynomial space. 

4. Interpolate f* from r/ F to £ F . 

After following these steps, each solution flux point will contain the nodal value of the optimal polynomial 
approximation of the upwind flux within the solution polynomial space. 

III.C. Residual Over-Integration 

The second over-integration method focuses on preventing aliasing errors by over-sampling the residual 
function. Using the interpolation and projection operators from section III. A, over-integration of the residual 
proceeds as follows: 

1. Interpolate the solution from to 

2. Use the base FR method at the nodes r/ 1 and r/ F to evaluate the residual function at r) l . 

3. Project the higher-order residual polynomial to the solution polynomial space. 

4. Interpolate the projected residual polynomial from r] 1 to £ : . 

5. Update the solution at using the de-aliased residual polynomial. 

In the previous algorithm, there is a redundancy proceeding through steps 4, 5, and restarting back at 
step 1. In step 4, the projected residual is interpolated from the quadrature points to the solution points for 
the sole purpose of evaluating the updated solution in step 5. Restarting back at step 1 just interpolates the 
updated solution right back to the quadrature points. We can instead skip the interpolations all together 
and perform the solution update at the quadrature points since this is simply a linear combination of the 
solution and residual polynomials. This removes any need for using the solution points since all operations 
are performed at the quadrature points. 

The key feature of this second over-integration method is that it is very easy to implement within an 
existing code since it only requires the addition of a single matrix multiplication to each time (sub-)step. 
The code is basically mislead into thinking it is running a simulation with the solution order set to the 
quadrature order. This way the code sets up the simulation with the solution points actually being the 
quadrature points and all supporting data structures are created accordingly. As the simulation is advanced 
in time, the residual function is evaluated as usual within each time step or Runge-Kutta stage. Just before 
the solution is updated after a time step or Runge-Kutta stage, the higher-order residual polynomial is 
projected down to the solution polynomial space. Finally, the solution is updated using this de-aliased 
residual polynomial before continuing to the next time step. 
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IV. Governing Equations 


The Euler equations describe the motion of an inviscid fluid such as the isentropic vortex that is the 
topic of this paper. The 2D Euler equations for a calorically perfect, compressible fluid can be written in 
differential form as 



P 


1 

c2 

8 

1 


pVy 

d 

PV X 

d 

A 

pvl 

d 

A 

PV x Vy 

dt 

1 

dx 

pV X Vy 

_v x (pE + p)_ 

dy 

PV 2 V 

Vy ( PE+ P )_ 


where p, ifo p , and E respectively refer to the density, velocity in coordinate direction *, pressure, and the 
total energy per unit mass of the system. Thermodynamic closure is found through the total energy equation 


pE = 


P 

7-1 


\p{ v l 


(46) 


where 7 is the ratio of specific heats for the fluid in question. Additionally, the ideal gas law 


P — pRgasE 


(47) 


is utilized to relate the temperature, T, to the pressure and density where R gas is the specific gas constant 
of the fluid. The speed of sound, a, is computed using the relation 


a = 



\Jl RgasT 


(48) 


In the present work, the fluid is assumed to be calorically perfect air where 7 = 1.4 and the dimensionalized 
specific gas constant is R gas = 287.15 J/(kg-K). 


V. Numerical Results 

In the discussion of the results we will frequently refer to the degree of the solution polynomial used for 
a simulation as well as the degree of the polynomial used for over-integration. Therefore, we need to define 
a few notational conventions that will be used in this section to refer to these quantities. Let Vn refer to 
a solution polynomial of degree n that nominally results in a (n + l) th -order method and let Qm refer to 
over-integrating using a polynomial of degree m. Combining these, we will use VnQm to refer to a Vn 
solution that is over-integrated to Qm. 

V.A. Isentropic Euler Vortex Problem 

The isentropic Euler vortex problem is commonly used. 16, 1 ' ’ 20, 36, 42—4 ' for testing the order of accuracy of 
a numerical method because it is easy to implement and the exact solution is known at all times. This is 
especially the case for high-order methods that are designed to produce minimal numerical dissipation. This 
problem is ideal for liigh-order methods because their theoretical accuracy should allow them to propagate 
the vortex with minimal entropy production for an indefinite amount of time. 

In Ref. 48, we made an attempt to create a definition for the isentropic Euler vortex problem that 
unifies many of the variations that exist in literature. While performing some / 1 -refinement studies with 
this problem, we encountered a few instabilities of varying strengths. We were able to determine that some 
of these instabilities were caused by poorly chosen boundary conditions and the size of the computational 
domain being too small. By zeroing the mean flow velocity to keep the vortex stationary, we were able 
to circumvent the problem with the boundary conditions and successfully complete an /i-refinement study 
showing the long term stability of the FR method. 

Whereas the /i-refinement study was successful, some of the lower-resolution simulations for the stationary 
vortex still showed signs of instabilities. Since this problem has been shown to exhibit aliasing errors that 
can affect the stability of the numerical simulation , 1 1,20,21 we identified these cases as good candidates for 
testing over-integration as a means of reducing aliasing errors and improving stability. For completeness, we 
will also explore the effects of over-integration on an extremely under-resolved stationary vortex as well as 
its general effects related to the nodal quadrature through an /i/p-refinement study. 
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Figure 3: Contours of the initial density profile for the stationary vortex problem. 


Additional instabilities can be introduced within the system by periodic boundaries and we need to isolate 
the effects of over-integration on possible aliasing instabilities. Keeping the vortex stationary is the only way 
to run the vortex problem for any significant amount of time without using periodic boundary conditions. 
Therefore, we will use the stationary version of the isentropic Euler vortex problem as defined in Ref. 48. 

A condensed definition for the stationary vortex problem follows since we are not concerned about unifying 
variations of the vortex problem found in literature. The initial conditions for the vortex problem are found 
by superimposing perturbations in velocity and temperature onto a uniform mean flow. These perturbations 
are defined such that the entropy, S = p/p 7 , is constant throughout the entire grid domain, i.e., 5 S = 0. 
In order to keep the vortex stationary, the mean flow velocity, (v Xt00 ,v yj00 ), is required to be zero and the 
free-stream density, p^, and pressure, p 0 OJ are set to one. 

The perturbations in velocity and temperature are defined as 


Sv x = - y Cl 
SVy = + X Cl 


(49) 


ST = — 




where Cl is a Gaussian function of the form 


n = A_ 

2tT\/7 


(50) 


These perturbations, along with the free-stream conditions and the isentropic relations between density, 
pressure, and temperature, are used to define the initial flow conditions of the primitive variables as 


Po = (1 + ST) 


Ar: . o — Sv x 


Vy^Q — 


(51) 


p 0 = A(i + «y 


where the subscript 0 denotes a quantity given at t = 0. The tilde accents on the primitive variables indicate 
that each is a nondimensional quantity where Poo, a^, and T, ' x have been used as the characteristic density, 
velocity, and temperature scales, respectively. 

The last remaining item needed to complete the problem is the definition of the computation do- 
main. For all of the simulations in this work, we used a square domain with grid boundaries located at 
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(x,y) = (±10, ±10). A number of different grid resolutions was used with each grid divided into regular, 
nonoverlapping quadrilateral cells of equal size. Finally, all grid boundaries were set to characteristic bound- 
ary conditions to allow for information to flow in/out of the domain as needed. A contour plot showing the 
initial density profile resulting from this problem definition is shown in figure 3. 

The beauty of this problem is that it can be run indefinitely and the exact analytical solution at any 
point in time is simply the initial conditions. For the propagating vortex, time can additionally be measured 
in flow-through periods or the number of times the vortex propagates completely through the domain and 
returns back to its initial location. Even though a flow-through period for the stationary vortex is undefined, 
we will still use it to measure the elapsed time of our simulations. In this case, one period is equivalent to 
the amount of time for a propagating vortex with v X}0a = 1 to complete one flow-through period. 

V.B. Over-Integration Applied to a Very Under- Resolved Case 

The first case used for testing over-integration is an attempt to visualize the effect that it has on the modes of 
the solution. Aliasing errors are more commonly found with under-resolved problems. In these cases, there 
is modal energy that should exist in modes higher than the span of the polynomial basis functions. These 
energies are erroneously resolved by the highest modes of the polynomial basis resulting in an aliased solution. 
Therefore, we should be able to use an extremely under-resolved case to see how well over-integration is able 
to remove energy from the highest modes in the areas where these modes should not be activated. 

To create a very under-resolved case, we used a grid containing 10 2 cells and a V8 solution with and 
without over-integration to V8QVd. The simulations were run for one period and the absolute value of the 
density modes were plotted for each grid cell in figure 4, with all modes less than 1.0 x 10 - ' cut off (set to 
white) . 

The two images at the top of figure 4 are given to illustrate how the modes are to be interpreted within 
each cell. These two images show up-close views of the grid cell bounded by the lines x = 0, x = 2, y = 2, 
and y = 4 from the respective image located directly below. The box in the lower left corner of the cell gives 
the value of the modal coefficient, c m , for the polynomial term c m x°y ° , i.e., the constant mode. Moving 
horizontally, the exponent for x increases from left to right. Moving vertically, the exponent for y increases 
from bottom to top. The box at the top right gives the value of the modal coefficient for the polynomial 
term c m x 8 y 8 . 

As discussed previously, we would expect to see the modes or boxes towards the top and right sides of 
each grid cell to be cut off or darker blue when using over-integration. More importantly, we would expect 
to see the boxes in the outer regions of the domain to have only the constant coefficient activated since 
those regions only contain the mean flow. Both of these expectations are clearly met when comparing the 
over-integrated modal solution at the bottom right of figure 4 to the base solution located at the bottom left 
of figure 4. This is a clear indication that over-integration is working as intended. 

V.C. Over-Integration Applied to Well-Resolved Cases 

We know that over-integration is working as intended in the previous under-resolved case, but what happens 
if we apply it to a moderately or well-resolved case? We tested various levels of over-integration on five 
different grid levels, 40 2 , 80 2 , 120 2 , 160 2 , and 200 2 , and using V2-V4 order solutions. For each of the 
combinations of grid level and polynomial order, we computed the density error through one period using 
three “base” methods without over-integration: Gauss-LP, Lobatto-LP, and Lobatto-CR. We then compare 
those results to simulations with varying degrees of over-integration applied to the Gauss-LP and Lobatto-LP 
methods. We did not use over-integration with the Lobatto-CR method since it already offers an additional 
order of accuracy over the Lobatto-LP method. Instead we wanted to see how over-integration with the 
Lobatto-LP method compares to the base Lobatto-LP and Lobatto-CR methods. 

There are a couple of consistencies between figure 5 and the following figures 6 and 7 to identify before 
continuing. First, across all three figures the dashed lines give the base Lobatto-LP (long dashes) and 
Lobatto-CR (short dashes) results and the solid line with square symbols gives the Gauss-LP results. Second, 
all the results on a given grid resolution are plotted using the same color in all the three figures. More 
specifically, the red, green, blue, magenta, and black groups of lines respectively represent the grid levels 
containing 40 2 , 80 2 , 120 2 , 160 2 , and 200 2 cells. 

Figure 5 shows the results for all of the V2 order simulations on each of the grid levels. Two sets of 
simulations involving over-integration are shown alongside the three base methods: Gauss-LP (solid line with 
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Figure 4: Comparison of the V8 modal solution without over-integration (left) and with over-integration to 
■P8Q16 (right). The top two images show a detailed view of the modes for the cells highlighted in red in the 
respective image below. All modes less than 1.0 x 10“ 7 have been cut off (set to white). 
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Period 

Figure 5: Comparison of V2 solutions with and without over-integration on multiple grid levels. 


‘x’ symbols) and Lobatto-LP (solid line with circle symbols) both over-integrating to V2Q4,. The results for 
the three base methods follow the same trends as found by other researchers: 33 the error of the Lobatto-LP 
method is between 0. 5-1.0 magnitude higher than the Gauss-LP method and the error of the Lobatto-CR 
method is just slightly higher than the Gauss-LP method. The very interesting result is that the error from 
both over-integration methods is identical to the base Gauss-LP method across all grid levels. 

Figure 6 shows the results for all of the 7^3 order simulations on each of the grid levels. As with the V2 
results, we have included two sets of simulations involving over-integration alongside the three base methods. 
Again, we have included over- integration to 7^3 Q6 with the Gauss-LP method (solid line with circle symbols) 
but this time we reduced the over-integration with the Lobatto-LP method (solid line with ‘x’ symbols) by 
one order to 7^325. Even with this change, all of the over-integration results are again identical to the base 
Gauss-LP method. 

Finally, figure 7 shows the results for all of the 7^4 order simulations on each of the grid levels. For 
this case, we did not include any over-integration results for the Gauss-LP method and instead focused on 
finding what degree of over-integration is required for the Lobatto-LP method to reach the same result as the 
base Gauss-LP method. We have included three sets of simulations with varying degrees of over-integration 
with the Lobatto-LP method alongside the three base methods: over-integrating to 7^428 (solid line with 
diamond symbols), 7^426 (solid line with circle symbols), and 7^425 (solid line with ‘x’ symbols). Again, all 
of the over-integration results are identical to the base Gauss-LP results. 

There are three important conclusions that can be made from applying over-integration to these mod- 
erately to well-resolved cases. First, when using Gauss solution points, over-integrating does nothing and 
gives the same result as if it had not been used at all. Second, when using Lobatto solution points, over- 
integration improves the answer to that found by using Gauss solution points. Finally, over-integrating by 
just one degree higher might be enough to remove any aliasing errors inherent to using Lobatto solution 
points. 
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Figure 6: Comparison of V3 solutions with and without over-integration on multiple grid levels 
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Figure 7: Comparison of P4 solutions with and without over- integration on multiple grid levels 
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V.D. Over-Integration Applied to Instabilities 

In Ref. 48, some under-resolved cases for the stationary vortex were found that contained instabilities of 
varying strengths. These instabilities can be split into two categories: “fatal” and “non-fatal” instabilities. 
A fatal instability is one that caused the simulation to crash due to negative pressure or density. Non-fatal 
instabilities are those cases that demonstrate instabilities but are still able to finish the simulation without 
crashing. What we want to see is whether over-integration can suppress these instabilities and provide some 
additional robustness when applied to an under-resolved simulation. We chose two cases for this test, one 
exhibiting a fatal instability and the second exhibiting a non-fatal instability. 

The case involving a fatal instability uses the Gauss-LP method with a V2 solution on a 40 2 grid and 
the base result without over-integration is seen in figure 8a as the black line with square symbols. We 
repeated this test case using the weak form of the DG Spectral Element Method (DGSEM) to ensure that 
this instability is not something that was only found using the FR method. The DGSEM results are shown 
as a blue line with diamond symbols and is exactly identical to the FR Gauss-LP results. Also shown in 
figure 8a are two additional results that use over-integration with the Gauss-LP method: the first over- 
integrating to V2QA (green line with ‘x’ symbols) and the second over-integrating to V2Q\0 (red line with 
circle symbols). In figure 8a, where each of the lines terminates is the point in time at which the simulation 
blows up. Without over-integration, the simulation blew up after four periods for both of the methods 
without over-integration. Adding over-integration was not able to improve the stability for this test case. 
Over- integrating to 7^2 Q4 slightly reduced the magnitude of the overall error but the simulation still diverged 
at the same point as the simulations without over- integration. Over- integrating to V2 Q10 reduced the error 
and allowed to simulation to proceed longer but only for an additional half period. 

The non-fatal instability case is the same as the previous case except that the grid resolution has been 
increased to a grid containing 80 2 cells. The base result using the Gauss-LP method is shown in figure 8b as 
the black line with square symbols. As done with the previous case, we over-integrated to V2Q1Q to ensure 
we had enough resolution to stabilize the solution. However, over-integration did not improve the stability 
of the simulation and gave a very similar solution to the simulation without over-integration. 

VI. Conclusions 

The results presented in this work give conflicting stories on the success of using over-integration as 
a means of preventing aliasing errors for the FR and DG methods. As shown in sections V.B and V.C, 
over-integrating indeed improves the quality of the solution, especially for under-resolved simulations and 
for methods utilizing a reduced order quadrature (e.g., solution points are Lobatto points). However, in 
section V.D, over-integrating applied to some test cases involving numerical instabilities does not improve 
the robustness of the FR and DG methods. 

Because of the difficulties that we have experienced with the vortex problem, 48 it seems likely that the 
instabilities we were trying to fix with over-integration in section V.D are not due to aliasing errors. One 
reason for this conjecture is a small detail from the final jump in error that is found in both simulations shown 
in figure 8b. These jumps coincide with an exponential fit found in Ref. 48 that appeared to be related to the 
accumulation of machine roundoff errors caused by simulating a closed system with all periodic boundaries. 
Conversely, this could also be a coincidence and there is still some sort of aliasing driven instability that has 
yet to be resolved. Either way, further testing is required to determine the cause of these instabilities that 
we have encountered with the stationary vortex problem. 

The results from section V.C show that over-integration successfully removes any aliasing errors that are 
inherent to the use of Lobatto quadratures. This is especially promising when looking to use triangular grid 
cells in 2D and tetrahedra and other unstructured cell types in 3D where Gaussian quadrature rules are 
not defined. Nodal sets based on a barycentric-type mapping of Lobatto quadrature for triangular elements 
is popular with many unstructured nodal DG methods. Our results indicate that over-integrating by an 
order or two might be sufficient in removing most of the aliasing errors found in the Lobatto nodal sets for 
triangles. In general when the flow is well-resolved, aliasing errors often do not cause numerical instabilities. 
The focus of future work includes applying over-integration to the Navier-Stokes equations and then to 3D 
unstructured cells (tetrahedra, pyramids, prisms, and hexahedra geometries). 
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P2 - 40x40 - Density Error 



(a) Over-integration applied to a fatal instability. Where each line terminates is the point that the simulation failed. 


P2 - 80x80 - Density Error 



(b) Over-integration applied to a non-fatal instability. The P2Q10 simulation terminates before it reaches 100 periods 
because it reached the total wall time the job had been allocated and not because it diverged. 

Figure 8: Application of different amounts of over-integration to reduce instabilities caused by insufficient 
resolution. 


20 of 22 


American Institute of Aeronautics and Astronautics 


References 


^^DeBonis, J., “Progress Towards Large- Eddy Simulations for Prediction of Realistic Nozzle Systems,” AIAA Journal of 
Propulsion and Power , Vol. 23, No. 5, 2007, pp. 971-980. 

2 Forsythe, J., Squires, K., Wurtzler, K., and Spalart, P., “Detached-Eddy Simulation of Fighter Aircraft at High-Alpha,” 
Journal of Aircraft, Vol. 41, No. 2, 2004, pp. 193-200. 

3 Mavripilis, D., Pelaez, J., and Kandil, O., “Large Eddy and Detached Eddy Simulations Using an Unstructured Multigrid 
Solver,” DNS/LES - Progress and Challenges. Proceedings of the Third AFOSR International Conference on DNS/LES , 
Columbus, OH, 2001. 

4 Moreau, S., Christophe, J., and Roger, M., “LES of the Trailing-Edge Flow and Noise of a NACA0012 Airfoil Near Stall,” 
Proceedings of the Summer Program 2008, Center for Turbulence Research , Stanford University/ NASA Ames, 2008. 

5 Kang, S., Iaccarino, G., Ham, F., and Moin, P., “Prediction of Wall-Pressure Fluctuation in Turbulent Flows with an 
Immersed Boundary Method,” Journal of Computational Physics , Vol. 228, No. 18, 2009, pp. 6753-6772. 

6 Huynh, H., “A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods,” AIAA 
Paper 2007-4079, Jun 2007. 

'Huynh, H., “A Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin for Diffusion,” AIAA 
Paper 2009-403, Jan 2009. 

8 Huynh, H., “High- Order Methods Including Discontinuous Galerkin by Reconstructions on Triangular Meshes,” AIAA 
Paper 2011-44, Jan 2011. 

9 Vermeire, B. C., Cagnone, J.-S., and Nadarajah, S., “ILES Using the Correction Procedure via Reconstruction Scheme,” 
AIAA Paper 2013-1001, Jan 2013. 

10 Vermeire, B. C., Nadarajah, S., and Tucker, P. G., “Canonical Test Cases for High-Order Unstructured Implicit Large 
Eddy Simulation,” AIAA Paper 2014-0935, Jan 2014. 

11 Skarolek, V. and Miyaji, K., “Transitional Flow over a SD7003 Wing Using Flux Reconstruction Scheme,” AIAA Paper 
2014-0250, Jan 2014. 

12 Haga, T., Tsutsumi, S., Kawai, S., and Takaki, R., “Large-Eddy Simulation of a Supersonic Jet Using High-Order Flux 
Reconstruction Scheme,” AIAA Paper 2015-0831, Jan 2015. 

13 Kirby, R. M. and Karniadakis, G. E., “De- Aliasing on Non-Uniform Grids: Algorithms and Applications,” Journal of 
Computational Physics , Vol. 191, No. 1, 2003, pp. 249—264. 

14 Kirby, R. M. and Sherwin, S. J., “Aliasing Errors Due to Quadratic Nonlinearities on Triangular Spectral/hp Element 
Discretisations,” Journal of Engineering Mathematics, Vol. 56, No. 3, 2006, pp. 273-288. 

15 Kopriva, D. A., “Metric Identities and the Discontinuous Spectral Element Method on Curvilinear Meshes,” Journal of 
Scientific Computing, Vol. 26, No. 3, 2006, pp. 301-327. 

16 Hesthaven, J. S. and Warburton, T., Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, 
Vol. 54 of Texts in Applied Mathematics, Springer New York, 2008. 

17 Castonguay, P., Vincent, P., and Jameson, A., “Application of High-Order Energy Stable Flux Reconstruction Schemes 
to the Euler Equations,” AIAA Paper 2013-686, Jan 2011. 

18 Gassner, G. J. and Beck, A. D., “On the Accuracy of High-Order Discretizations for Underresolved Turbulence Simula- 
tions,” Theoretical and Computational Fluid Dynamics, Vol. 27, No. 3-4, 2013, pp. 221-237. 

19 Karniadakis, G. and Sherwin, S., Spectral/hp Element Methods for Computational Fluid Dynamics, Oxford University 
Press, 2nd ed., 2013. 

20 Williams, D. M. and Jameson, A., “Nodal Points and the Nonlinear Stability of High-Order Methods for Unsteady Flow 
Problems on Tetrahedral Meshes,” AIAA Paper 2013-2830, Jun 2013. 

21 Witherden, F. D. and Vincent, P. E., “An Analysis of Solution Point Coordinates for Flux Reconstruction Schemes on 
Triangular Elements,” Journal of Scientific Computing, Vol. TBD, No. TBD, 2014, pp. TBD. 

22 Kanevsky, A., Carpenter, M. H., and Hesthaven, J. S., “Idempotent Filtering in Spectral and Spectral Element Methods,” 
Journal of Computational Physics, Vol. 220, No. 1, 2006, pp. 41-58. 

23 Tadmor, E., “Convergence of Spectral Methods for Nonlinear Conservation Laws,” SIAM Journal on Numerical Analysis, 
Vol. 26, No. 1, 1989, pp. 30-44. 

24 Maday, Y., Ould Kaber, S. M., and Tadmor, E., “Legendre Pseudospectral Viscosity Method for Nonlinear Conservation 
Laws,” SIAM Journal on Numerical Analysis, Vol. 30, No. 2, 1993, pp. 321-342. 

25 Karamanos, G.-S. and Karniadakis, G. E., “A Spectral Vanishing Viscosity Method for Large-Eddy Simulations,” Journal 
of Computational Physics, Vol. 163, No. 1, 2000, pp. 22-50. 

26 Guo, B.-y., Ma, H.-p., and Tadmor, E., “Spectral Vanishing Viscosity Method for Nonlinear Conservation Laws,” SIAM 
Journal on Numerical Analysis, Vol. 39, No. 4, 2001, pp. 1254-1268. 

27 Kirby, R. M., Toward Dynamic Spectral/hp Refinement: Algorithms and Applications to Flow- Structure Interactions, 
Ph.D. thesis, Brown University, Providence, RI, 2003. 

28 Funaro, D., “Some Remarks about the Collocation Method on a Modified Legendre Grid,” Computers & Mathematics 
with Applications, Vol. 33, No. 1, 1997, pp. 95-103. 

29 Funaro, D., Spectral Elements for Transport- Dominated Equations, Springer Berlin Heidelberg, 1997. 

30 Vincent, P., Castonguay, P., and Jameson, A., “A New Class of High-Order Energy Stable Flux Reconstruction Schemes,” 
Journal of Scientific Computing, Vol. 47, No. 1, 2011, pp. 50-72. 

31 Jameson, A., Castonguay, P., and Vincent, P., “On the Non-Linear Stability of Flux Reconstruction Schemes,” Journal 
of Scientific Computing, Vol. 50, No. 2, 2012, pp. 434-445. 

32 De Grazia, D., Mengaldo, G., Moxey, D., Vincent, P. E., and Sherwin, S. J., “Connections Between the Discontinuous 
Galerkin Method and High-Order Flux Reconstruction Schemes,” International Journal for Numerical Methods in Fluids , 
Vol. 75, No. 12, 2014, pp. 860-877. 


21 of 22 


American Institute of Aeronautics and Astronautics 


33 Yu, M., Wang, Z., and Liu, Y., “On the Accuracy and Efficiency of Discontinuous Galerkin, Spectral Difference and 
Correction Procedure via Reconstruction Methods,” Journal of Computational Physics , Vol. 259, 2014, pp. 70-95. 

34 Kopriva, D. A., Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists and Engi- 
neers , Springer Science & Business Media, 2009. 

35 Wang, Z. and Gao, H., “A Unifying Lifting Collocation Penalty Formulation Including the Discontinuous Galerkin, 
Spectral Volume/ Difference Methods for Conservation Laws on Mixed Grids,” Journal of Computational Physics , Vol. 228, 
No. 21, 2009, pp. 8161-8186. 

36 Gao, H. and Wang, Z., “A Conservative Correction Procedure via Reconstruction Formulation with the Chain-Rule 
Divergence Evaluation,” J. Comput. Physics , Vol. 232, No. 1, 2013, pp. 7-13. 

37 Roe, P. L., “Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes,” Journal of Computational 
Physics , Vol. 43, No. 2, 1981, pp. 357-372. 

38 Huynh, H. T., “Accurate Upwind Methods for the Euler Equations,” SIAM Journal on Numerical Analysis , Vol. 32, 
No. 5, 1995, pp. 1565-1619. 

39 Gottlieb, S. and Shu, C.-W., “Total Variation Diminishing Runge-Kutta Schemes,” Mathematics of Computation , Vol. 67, 
No. 221, 1998, pp. 73-85. 

40 Kopriva, D. A. and Gassner, G., “On the Quadrature and Weak Form Choices in Collocation Type Discontinuous Galerkin 
Spectral Element Methods,” Journal of Scientific Computing , Vol. 44, No. 2, 2010, pp. 136-155. 

41 Gassner, G. and Kopriva, D. A., “A Comparison of the Dispersion and Dissipation Errors of Gauss and Gauss-Lobatto 
Discontinuous Galerkin Spectral Element Methods,” SIAM Journal on Scientific Computing , Vol. 33, No. 5, 2011, pp. 2560- 
2579. 

42 Shu, C.-W., “Essentially Non-oscillatory and Weighted Essentially Non-oscillatory Schemes for Hyperbolic Conservation 
Laws,” Advanced Numerical Approximation of Nonlinear Hyperbolic Equations , edited by A. Quarteroni, Vol. 1697 of Lecture 
Notes in Mathematics , Springer Berlin Heidelberg, 1998, pp. 325-432. 

43 Wang, Z. J., Liu, Y., May, G., and Jameson, A., “Spectral Difference Method for Unstructured Grids II: Extension to 
the Euler Equations,” Journal of Scientific Computing , Vol. 32, No. 1, 2007, pp. 45-71. 

44 Vincent, P., Castonguay, P., and Jameson, A., “Insights from von Neumann Analysis of High-Order Flux Reconstruction 
Schemes,” Journal of Computational Physics , Vol. 230, No. 22, 2011, pp. 8134-8154. 

45 Wang, Z., Fidkowski, K., Abgrall, R., Bassi, F., Caraeni, D., Cary, A., Deconinck, H., Hartmann, R., Hillewaert, K., 
Huynh, H., Kroll, N., May, G., Persson, P.-O., van Leer, B., and Visbal, M., “High-order CFD Methods: Current Status and 
Perspective,” International Journal for Numerical Methods in Fluids , Vol. 72, No. 8, 2013, pp. 811-845. 

46 Witherden, F. D., Farrington, A. M., and Vincent, P. E., “PyFR: An Open Source Framework for Solving Advection- 
Diffusion Type Problems on Streaming Architectures Using the Flux Reconstruction Approach,” Computer Physics Commu- 
nications , Vol. 185, No. 11, 2014, pp. 3028-3040. 

47 Yu, M. and Wang, Z. J., “On the Accuracy and Efficiency of Several Discontinuous High-Order Formulations,” AIAA 
Paper 2013-855, Jan 2013. 

48 Spiegel, S., Huynh, H. T., and DeBonis, J. R., “Survey of the Isentropic Euler Vortex Problem Using High-Order 
Methods,” AIAA Paper 2015-TBD, Jun 2015. 


22 of 22 


American Institute of Aeronautics and Astronautics 


